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ABSTRACT 


Small  underground  nuclear  explosions  need  to  be  confidently  detected  and  identified  in  regions  of  the  world  where 
they  have  never  occurred.  We  are  developing  a  parametric  model  of  the  nuclear  explosion  seismic  source  spectrum 
derived  from  regional  phases  ( Pn ,  Pg,  Sn,  and  Lg )  that  is  compatible  with  earthquake-based  geometrical  spreading 
and  attenuation.  Earthquake  spectra  are  fit  with  a  generalized  version  of  the  Brune  spectrum,  which  is  a  three- 
parameter  model  that  describes  the  long-period  level,  corner-frequency,  and  spectral  slope  at  high  frequencies. 

These  parameters  are  then  correlated  with  near-source  geology  and  containment  conditions.  There  is  a  correlation  of 
high  gas-porosity  (low  strength)  with  increased  spectral  slope.  However,  there  are  trade-offs  between  the  slope  and 
comer- frequency,  which  we  try  to  independently  constrain  using  Mueller-Murphy  relations  and 
coda-ratio  techniques.  The  relationship  between  the  parametric  equations  and  the  geologic  and  containment 
conditions  will  assist  in  our  physical  understanding  of  the  nuclear  explosion  source.  The  achievable  goal  of  our 
parametric  model  development  is  to  be  able  to  predict  observed  local  and  regional  distance  seismic  amplitudes  for 
event  identification  and  yield  determination  in  regions  with  incomplete  or  no  prior  history  of  underground  nuclear 
testing. 
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OBJECTIVES 


We  aim  to  develop  a  practical  explosion  source  parametric  spectral  model,  based  on  all  available  data,  that  describes 
nuclear  explosion  P-  and  S-wave  source  spectra  for  a  variety  of  geologic  and  containment  conditions.  This  approach 
follows  the  simple  earthquake  parametric  spectral  model  based  on  Brune  (1970),  which  is  used  for  the  MDAC 
approach  (Walter  and  Taylor,  2002)  to  improve  earthquake/explosion  discrimination.  In  regions  without  prior 
explosions,  the  parametric  model  could  be  combined  with  earthquake- derived  path  corrections  to  predict  explosion 
regional  phase  amplitudes,  improve  discriminants  such  as  P/S  ratios,  and  support  identification  procedures 
(e.g.,  Event  Classification  Matrix,  [ECM])  that  explicitly  need  to  use  explosion  discriminant  probability  density 
functions. 

It  is  well  known  that  depth  and  near-source  material  properties  can  affect  seismic  estimates  of  explosion  yield,  and 
prior  work  at  the  Nevada  Test  Site  (NTS)  (e.g.,  Walter  et  al.,  1995)  has  found  that  explosions  in  weak  materials  have 
lower  corner  frequencies  and  steeper  spectral  fall-offs  for  P- waves  than  is  predicted  by  the  standard  Mueller  and 
Murphy  (1971)  model.  As  part  of  this  research,  we  hope  to  quantify  these  effects  as  a  function  of  frequency  and 
wave  type.  Additionally,  many  of  the  most  effective  regional  discriminants  (high-frequency  P/S  ratios)  make  use  of 
S- waves,  as  do  S'- wave  coda  yield  estimation  techniques,  yet  there  remain  many  questions  about  how  to  predict 
explosion  S-wave  amplitudes.  The  development  of  a  combined  P-  and  S-wave  spectral  model  consistent  with 
observed  regional  P- and  S-wave  data  is  a  goal  of  this  work. 

RESEARCH  ACCOMPLISHED 


Introduction 


In  previous  work  that  looked  at  low  to  high  frequency  ratios  of  regional  phase  (e.g.,  Pn ,  Pg ,  Lg)  amplitudes  to 
separate  explosions  from  earthquakes  at  NTS,  Walter  et  al.  (1995)  noted  that  the  results  showed  a  strong  dependence 
on  the  source  media  properties.  Nuclear  tests  in  weak  and/or  high  gas  porosity  media  tended  to  have  higher  values 
and  discriminate  better  from  earthquakes  than  explosions  in  stronger  and/or  lower  gas  porosity  media. 


a)  MNV-KNB  average  Lg  spectral  ratios  and  model  fits 


b)  Source  spectrum 


c  )  Lg  transfer  function 


Ml  (coda) 


Figure  1.  Simple  first-order  source  model  fits  to  low  to  high  frequency  spectral  ratios  of  Lg  amplitudes  from 
Walter  et  al.  (1995,  Figure  8).  Earthquake  displacement  spectra  are  fit  with  the  Brune  (1970)  model 
(blue  line),  which  is  constant  at  low  frequencies  and  falls  off  above  a  corner  frequency  as /'2. 
Explosions  are  fit  with  two  extremes  of  the  Denny  and  Goodman  (1990)  model  in  order  to 
investigate  explosion  dependence  on  emplacement  material  (namely,  gas  porosity  [Gp]).  The  ‘Low 
Gp 9  model  (black  line)  has  an  effective  fall-off  similar  to  the  earthquake  model.  The  ‘High  Gp 9 
model  (red  line)  has  a  greater  fall-off  of /~3. 
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Figure  2.  a)  Regional  map  of  stations  used  in  the 

analysis.  The  location  of  the  map  within  the 
continental  US  is  given  by  the  inset  map. 
The  NTS  is  outlined  and  the  shaded  section 
is  shown  in  (b).  b)  Map  of  the  northern 
NTS  with  explosion  locations  and  Vergino 
and  Mensing  (1990)  area  designations. 


To  fit  the  explosion  data,  we  used  two  extremes  of  the 
Denny  and  Goodman  explosion  model  (1990),  which  has 
two  corner  frequencies.  In  the  low  Gp  case  we  allowed 
the  second  comer  to  be  at  a  higher  frequency  than  the 
range  of  interest,  giving  an  effective /"2  falloff.  In  the  high 
Gp  case  we  forced  the  comers  to  be  the  same,  giving  an 
/"3  falloff.  In  both  cases  we  used  the  observed  3  Hz  comer 
frequency  of  the  1993  Non-Proliferation  Experiment 
(NPE),  a  kiloton  chemical  explosion  and  assumed  the 
comer  frequency  scales  with  the  cube  root  of  ML( coda). 
Given  that  a  pure  explosion  should  not  generate  S-waves, 
one  way  to  think  about  the  Lg  spectral  ratio  is  as  the 
product  of  the  P-wave  source  ratio  and  a  transfer  function 
ratio,  where  the  transfer  function  is  a  representation  of 
how  efficiently  the  source  generated  P- waves  are 
converted  (by  whatever  means)  into  S-waves.  We 
estimated  the  transfer  function  ratio  as  a  function  of 
frequency  as  shown  in  Figure  1 ,  and  then  multiplied  the 
P- wave  based  Denny-Goodman  model  curves  by  these 
factors  and  then  compared  them  to  the  Lg  spectral  ratio  in 
Figure  1.  The  result  is  a  fairly  reasonable  first  order  fit  to 


the  data.  Interestingly,  since  the  explosion  Lg  transfer 

functions  are  larger  than  those  for  the  earthquakes,  it  shifts  the  explosion  Lg  spectral  ratios  to  relatively  higher 
values  than  for  the  earthquakes,  improving  the  discrimination  performance  of  Lg  spectral  ratios  over  those  of  Pn 
alone  (Walter  et  al.  1995,  Figure  7). 


For  example,  in  Figure  1  the  ratio  of  the  Lg  amplitude  at 
1-2  Hz  compared  with  the  amplitude  at  6-8  Hz  is  shown 
as  a  function  of  magnitude.  The  earthquakes  (blue  circles 
and  green  squares)  show  the  expected  trend  with 
magnitude  going  from  a  low  value  for  small  events  when 
the  comer  frequency  is  above  8  Hz  and  both 
measurements  are  on  the  constant  part  of  the  source 
spectra  that  is  proportional  to  moment.  For  large 
magnitudes,  the  source  comer  frequency  drops  below  1 
Hz  and  then  both  measures  are  on  the  part  of  the  source 
spectra  that  decays  with  frequency  as /"2,  resulting  in  high 
spectral  ratio  values.  As  magnitude  increases,  the 
earthquakes  follow  a  sigmoid  curve  as  shown  by  the  blue 
line,  which  is  based  on  the  Bmne  (1970)  model.  The 
explosions  are  split  into  two  categories  based  on  the 
source  media,  a  high  gas  porosity  ( Gp )  and  low  strength 
(defined  as  pot,  where  p  is  density  and  a  is  the  local 
compressional  wave  speed)  group  (red  x)  and  low  Gp 
high  strength  group  (black  crosses),  and  a  clear  difference 
between  the  two  can  be  seen.  In  fact,  the  low  Gp 
explosions  reach  spectral  ratio  values  that  imply  much 
steeper  falloff  than  f2. 


While  in  this  case,  the  Denny  and  Goodman  (1990)  model 
provides  a  reasonable  fit  to  these  NTS  data,  it  is  not  clear  how  we 
would  make  use  of  it  in  other  regions.  Furthermore  the  limitation 
of  only  two  choices  of  high  frequency  fall-off  does  not  capture  the 
full  range  of  the  observations.  To  develop  a  practical  parametric 
source  spectral  model  we  want  to  be  able  to  tie  parameters  like 
low-frequency  level,  comer- frequency  and  fall-off  rate  to 


Table  1.  NTS  regional  attenuation  and 
_ geometric  spreading  parameters. 

Phase  p  r0  (km)  Q0  y 


Pn 

1.10.  001 

210  0. 

65 

pg 

0.5  1  00 

190  0 

.45 

Lg 

0.5  1  00 

200  0 

.54 
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measureable  properties  like  yield,  depth,  and  media  properties 
such  as  gas  porosity,  water  content,  and  strength.  In  the  next 
section  we  take  a  more  general  approach  to  fitting  the  NTS 
explosion  data  with  a  simple  spectral  model  parameterized  by 
a  long-period  level,  comer- frequency,  and  high-frequency  roll¬ 
off.  We  then  relate  these  parameters  for  Pn  spectra  to  source 
and  medium  properties.  In  the  future  we  will  perform  similar 
analyses  for  Pg  and  Lg  spectra. 

Data  and  Methods 

We  employ  the  NTS  explosion  dataset  of  Walter  et  al.  (2004), 
specifically  the  raw  spectra  of  that  dataset.  Waveforms  are 
de-meaned,  de-trended  and  instmment  corrected  to 
acceleration.  The  signal  is  windowed  with  a  5%  cosine  taper 
that  starts  before  the  pick,  where  the  time  before  the  pick  is  5% 
of  the  total  time  window.  The  Fourier  transform  is  calculated 
and  displacement  spectra  are  obtained  via  double  integration  in 
the  frequency  domain.  Finally,  the  resultant  amplitude  spectra 
are  interpolated  and  smoothed  to  obtain  a  sampling  period  of 
0.05  logio  Hz. 


Spectra  of  each  seismic  phase  for  each  explosion  are  calculated 
from  the  recordings  of  stations  of  the  Livermore  NTS 
Network,  ELK  (Elko,  NV),  KNB  (Kanab,  UT),  LAC  (Landers, 
CA),  and  MNV  (Mina,  NV).  The  locations  of  these  stations 
relative  to  the  NTS  are  given  in  Figure  2.  These  spectra  are 
then  corrected  for  geometrical  spreading  and  regional, 
frequency-dependent  attenuation  of  the  form  Q  =  Qofr,  where 
Qo  is  Q  at  1  Hz  and  yis  the  power-law  dependence  on 
frequency,/  We  employ  the  Street  et  al.  (1975)  parametric 
form  of  geometrical  spreading, 


G(r)  = 


r  <ro 

r>r0 


(1) 


where  r0  is  the  distance  at  which  the  spreading  transitions  from 
spherical-  to  a  cylindrical-type  spreading  and  77  is  the  distance 
dependence.  The  attenuation  and  spreading  model  parameters 
for  each  seismic  phase  are  given  in  Table  1. 

We  require  three  of  the  four  stations  listed  above  to  have 
recorded  an  event  with  a  signal-to-noise  ratio  greater  than  two. 
We  also  inspect  each  spectra  for  non-stationary  noise  at  high- 
frequency  (typically  >9  Hz)  due  to  multi-band  recording 
problems  and  restrict  spectral  fitting  to  bands  unaffected  by 
problems.  The  spectra  are  jointly  fit  with  a  simple  parametric 
form  given  by 


S(f)-- 


1 +C///J 


(2) 


in  a  least-squares  inversion  that  also  provides  standard  error 
for  each  parameter  estimate.  Equation  (2)  describes  the 


Pn  spectra 


ie-03 


ic-06 


Ic-07 


le-0* 


le-07 


lc-QH 


!e~09 


Frequency  (Hz) 


Figure  3.  Pn  spectra  examples  with  the  name, 
magnitude  (where  M  =  mb[Pw]), 
depth  (Z),  and  gas  porosity  ( Gp ) 
given  in  the  subtitles.  The  long- 
period  level  (£#),  corner-frequency 
(fc),  and  fall-off  (y/)  are  given  in  each 
plot,  and  the  spectra  are  colored  by 
station  and  the  fit  is  shown  in  the 
legend  in  the  lowermost  plot, 
a)  Spectra  for  an  event  with  low  Gp 
and  co2  fall-off.  b)  Spectra  for  an 
event  with  high  Gp  and  high  fall-off. 
c)  Spectra  for  an  event  with  low  Gp 
and  high  fall-off. 
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Figure  4.  Regional  map  of  station  in  Borovoye  (BRVK)  and 
recorded  events  in  the  archive  (dots)  along  with 
the  event  shown  in  Figure  5  (star).  The  location  of 
the  map  is  given  in  the  inset  global  map. 

simplest  behavior  expected  for  seismic  spectra.  It  has  a  constant 
level  at  low  frequencies  (S0),  which  is  proportional  to  static 
displacement,  and  falls  off  at  high  frequencies  with  a  slope  of f~¥ 
beyond  a  comer  frequency  fc.  Examples  of  the  spectra  and  model 
fits  are  proportional  to  static  displacement,  and  fall  off  at  high 
frequencies  with  a  slope  of f~ ^  beyond  a  comer  frequency  fc. 
Examples  of  the  spectra  and  model  fits  are  given  in  Figure  3. 
BULLION  (Figure  3a)  is  fit  well  with  a  standard f'2  (y/=  2) 
spectral  fall-off,  but  MONTEREY  (Figure  3b),  detonated  in  weak 
material,  requires  a  steeper  fall-off  where  the  best  fit  y  «  4. 

We  plan  to  complement  the  NTS  dataset  with  the  Borovoye  (BRV) 
archive  that  has  recently  been  deglitched  and  response  functions 
estimated  (Richards  and  Kim,  2009).  This  archive  could  potentially 
provide  200  recordings  of  explosions  at  Semipalatinsk  Test  Site  of 
the  former  Soviet  Union  (Figure  4).  As  an  example,  we  plot  the 
regional  phase  spectra  of  one  of  the  explosions  mapped  in  Figure  4 
(shown  by  a  star)  in  Figure  5.  The  attenuation  and  spreading  model 
parameters  used  to  correct  these  spectra  are  given  in  Table  2.  The 
Pn  and  Pg  spectra  are  fit  with  y  «  3,  whereas  Lg  fall-off  is 
approximately  y  «  2. 


The  BRV  regional  seismic  phase  spectra,  especially  Pn  and  Pg 
(Figure  5a-b),  show  the  trade-off  in  fitting  the  comer  frequency 
and  the  roll-off.  In  this  case  a  slightly  smaller  fc  would  lead 
perhaps  a  more  appropriate  y  of  2.  Attempts  to  diminish  this  trade¬ 
off  will  be  discussed  in  the  next  section. 

Preliminary  results 

We  present  the  results  of  fitting  the  Pn  spectra  for  NTS  explosions 
in  Figure  6,  where  corner- frequency  is  plotted  versus  fall-off.  There 
is  a  correlation  between  the  two  estimated  parameters  that  is  most 
probably  related  to  the  least-squares  fitting  of  Equation  (2)  and  not 
to  any  physical  relationship.  This  trade-off  can  be  seen  graphically 


BRVK  spectra  (4  Oct  89  M4.7) 


Figure  5.  Pn  spectra  of  example  BRVK 

event  (shown  by  star  in  Figure  4). 
Values  given  in  each  plot  are 
defined  in  Figure  3. 


Table  2.  BRV  regional  attenuation  and 
_ geometric  spreading  parameters 

Phase  ij  r0  (km)  g0  y 


Pn 

1.10.  001 

210  0. 

65 

pg 

0.5  1  00 

190  0 

.45 

Lg 

0.5  1  00 

200  0 

.54 
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Pn  Pn 


Figure  6.  Correlation  of  estimated  source  parameters, 
corner-frequency  (fc )  and  fall-off  ( i//)  from 
Equation  (2). 

in  the  BRV  spectra  plotted  in  Figure  5.  The  trade-off 
could  produce  biased  estimates  of  y/,  so  we  seek  to 
constrain  fc  in  Equation  (2)  using  the  comer  frequency 
predicted  by  Mueller  and  Murphy  (1971)  as 
parameterized  by  Stevens  and  Day  (1985)  for 
tuff/rhyolite.  Figure  7a  shows  that  the  freely-fit  fc  is  well- 
correlated  with  the  Mueller  and  Murphy  (1971)  predicted 
comer- frequency  (MM  fc)  and  the  slope  is  close  to  one.  It 
can  be  made  even  closer  by  decreasing  the  reference 
elastic  radius  used  to  calculate  the  proportionality 
constant  used  in  the  Mueller  and  Murphy  (1971) 
relationship.  We  set/,  =  MM  fc  in  Equation  (2)  and 
estimate  y/.  This  new  y/is  less  correlated  with  MM  fc  as 
shown  in  Figure  7b. 

A  least-squares  regression  of  spectral  fall-off  ( y/)  for  Gp 
is  given  in  Figure  8a,  where  the  weights  are  given  by  the 
inverse  variance  of  ^/estimated  by  the  inversion.  Figure 
8b  is  the  same,  but  uses  the  estimated  y/  while  fixing  fc  as 
described  previously.  We  also  perform  a  comprehensive 


Equation  (2)  compared  with  those 
predicted  by  Mueller  and  Murphy  (1971), 
MM71  fc.  The  linear  relationship  is  given 
by  the  equation  at  the  top  of  the  plot  and 
the  slope  is  close  to  one.  b)  ^obtained 
while  fixing  fc  in  Equation  (2)  to  the  MM71 
fc.  The  correlation  of  the  terms  is  reduced 
compared  with  Figure  6. 


search  for  material  property  correlation  with  y/,  using  several  regressions  on  ^for  water  content  (Sw),  gas  porosity 
(Gp),  density  ( p ),  compressional  wave  speed  (a),  compressional  modulus  (pa),  strength  (pc?),  depth,  depth  to 
water-table,  and  scaled  depth-of-burial  (sDOB  =  depth/yield173).  These  geophysical  parameters  are  taken  from  the 
database  of  Springer  et  al.  (2002)  and  are  reported  to  represent  an  integrated  value  over  the  source  region.  Yield  was 
inferred  from  mb(Pn)  based  on  the  Vergino  and  Mensing  (1990)  relationship  at  NTS.  Using  a  step-wise  regression 
criteria  we  found  the  most  significant  parameters  in  a  linear  regression  model  to  be  Sw,  Gp,  and  sDOB. 

We  then  perform  a  least-squares  regression  using  these  parameters  weighted  by  the  estimated  inverse  variance  of  y/. 
The  model  is  given  by 

y/=  2.977  -4.026Gp+  3.5 19SW  +  2.81  x  10  ^  sDOB .  (3) 
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Figure  8.  Fall-off  (i//)  as  a  function  of  gas  porosity  (Gp)  for  a)  unconstrained  fc  and,  b)  fc  fixed  to  the  MM71- 
predicted  value.  The  best-fit  relationship  weighted  by  the  inverse  variance  of  the  estimated  ^is 
given  in  the  inset  equation. 

When  we  regress  using  ^estimated  while  fixing  fc,  we  get  a  slightly  different  model  given  by 

^  =  1.823  +  0.433 Gp  +  1.710^  +  2.2 1  x  10  3  sDOB.  (4) 

In  new  testing  environments  geophysical  parameters  such  as  gas  porosity  and  water  content  may  not  be  accessible. 
Using  NTS  geophysical  data  from  Springer  et  al.  (2002)  we  find  that  Gp  is  best  approximated  by  p ,  a ,  and  their 
interaction,  pa,  also  known  as  the  compressional  modulus.  This  relationship  is  given  by 

Gp  =  1 .266  -  5.437  x  10 ~4 p- 3.716  x  lO^a  +  l .635  x  10 ~7  pa,  (5) 

where  is  p  is  in  kg/m3  and  a  is  in  m/s. 

CONCLUSIONS  AND  RECOMMENDATIONS 

Our  preliminary  analysis  of  Pn  spectra  of  NTS  explosions  shows  that  a  simple  spectral  model  with  variable  fall-off 
is  appropriate,  and  the  fall-off  correlates  well  with  material  properties,  gas  porosity  and  water-content,  and  the 
scaled  depth-of-burial.  Comer- frequencies  obtained  with  the  variable  fall-off  model  correlate  with  Mueller  and 
Murphy  (1971)  corner- frequencies.  If  we  use  a  smaller  elastic  radius  than  the  one  given  for  tuff,  we  believe  we 
could  match  the  trend  of  the  Mueller-Murphy  comer  frequencies  with  a  slope  close  to  one. 
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Future  Work 

This  initial  P-wave  spectral  fitting  focused  on  correlating  material  properties  with  the  model  parameters.  We  will 
perform  similar  analyses  for  NTS  seismic  phases,  Pg  and  Lg ,  and  compare  the  results  with  those  presented  above  for 
Pn.  We  will  make  similar  calculations  for  spectra  at  Borovoye.  The  comparison  of  seismic  phases  will  allow  for  an 
examination  of  P/S  scaling  relationships,  which  can  be  compared  with  results  from  the  Fisk  (2007)  study  at  NTS. 
The  results  of  this  study  with  a  variable  fall-off  spectral  model  will  be  compared  with  other  spectral  models 
(e.g.,  Denny  and  Johnson  [1991]  and  references  therein).  In  addition,  we  will  formally  analyze  the  error  and  trade¬ 
offs  in  the  estimated  parameters  of  Equation  (2). 
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